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We numerically study a simple fluid composed of particles having a hard-core repulsion, comple- 
mented by two short-ranged attractive (sticky) spots at the particle poles, which provides a simple 
model for equilibrium polymerization of linear chains. The simplicity of the model allows for a close 
comparison, with no fitting parameters, between simulations and theoretical predictions based on 
the Wertheim perturbation theory, a unique framework for the analytic prediction of the proper- 
ties of self-assembling particle systems in terms of molecular parameter and liquid state correlation 
functions. This theory has not been subjected to stringent tests against simulation data for ordering 
across the polymerization transition. We numerically determine many of the thermodynamic prop- 
erties governing this basic form of self-assembly (energy per particle, order parameter or average 
fraction of particles in the associated state, average chain length, chain length distribution, average 
end-to-end distance of the chains, and the static structure factor) and find that predictions of the 
Wertheim theory accord remarkably well with the simulation results. 
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I. INTRODUCTION 

Recently, there has been great interest in exploiting 
self-assembly to create functional nanostructures in man- 
ufacturing, and this challenge has stimulated a great 
deal of experimental and theoretical activity [l], 0, 0, 
1, i, |, 0, B 0, M El El, El, El- Self-assembly has 
been considered for over 50 years to be central to un- 
derstanding structure formation in living systems and 
modeling and measurements of naturally occurring self- 
assembling systems has long been pursued in the biolog- 
ical sciences [Hj]. Even the term self-assembly derives 
from an appreciation of the capacity of viruses to spon- 
taneously reconstitute themselves from their molecular 
components [TBI [Trij , much as in the familiar example of 
micelle formation by block copolymers, lipids, and other 
surfactant molecules exhibiting amphiphilic interactions. 
The diversity and morphological and functional complex- 
ity of viruses, and the vast number of biological processes 
that form by a similar process in living systems, point to 
the potential of this type organizational process for man- 
ufacturing new materials. While the potential of self- 
assembly as a manufacturing strategy is clear, our un- 
derstanding of how this process actually works is still in- 
complete and many of the basic principles governing this 
type of organization are unclear. An evolutionary (trial 
and error) approach to this problem is not very efficient 
for manufacturing. There is evidently a need for develop- 
ing a first principles understanding of this phenomenon 
where no free parameters are involved in the theoreti- 
cal description to elucidate the fundamental principles 
governing self-assembly and the observables required to 
characterize the interactions governing thermodynamic 
self-assembly transitions, at least for simple model sys- 
tems that can be subjected to high resolution investiga- 



tion. 

As a starting point for this type of fundamental inves- 
tigation of self-assembly, we consider the problem of the 
equilibrium polymerization of linear polymer chains fl7l . 
EH, Ell, which is arguably the simplist variety of ther- 
mally reversible particle assembly into extended objects 
(polymers in the case of our molecular model). To inves- 
tigate this problem analytically, we exploit the Wertheim 
thermodynamic perturbation theory (W-TPT), which of- 
fers a systematic molecularly-based framework for calcu- 
lating the thermodynamic properties of self-assembling 
systems, although this theory has rarely before been ap- 
plied to this purpose [20( . The Wertheim theory has been 
previously considered to better understand properties of 
associating fluids and a similar sticky sphere model to 
the one we consider below has been considered to deter- 
mine how particle association affects critical properties 
associated with fluid phase separation (critical temper- 
ature and composition, binodals, critical compressibility 
factor, etc. [2l|, [22|, |23],_to determine the effect of asso- 
ciation on nucleation [24j and model antigen-antibody 
bonding [25| . In contrast, we are concerned here with 
the thermodynamic transition that accompanies the self- 
assembly of the particles into organized structures due to 
their anisotropic interactions. 

The Wertheim theory is certainly not a unique theory 
of the thermodynamics of self- assembling particle sys- 
tems. Models of equilibrium polymerization of linear, 
branched and compact structures have all been intro- 
duced based on the concept of an association equilibrium 
being established between the assembling particles and 
comparison of this type of theor y to simulation has led 
to remarkably good agreement |i~3l I2H |27l . l28l l29l . l3(il . 
HH, H3, Hl| • Up to the present time, however, it has been 
necessary to adjust the entropy of association parameter 
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FIG. 1: Pictorial representation of the model studied. Par- 
ticles are modeled as hard-core spheres (grey large sphere of 
diameter a), decorated by two sites located on the surface 
along a diameter (centers of the small gold spheres of diame- 
ter 25). The volume of the gold sphere outside the grey large 
sphere is the bonding volume. When the site of another parti- 
cles is located within the bonding volume, the pair interaction 
energy is taken equal to — Uo. 



units of the potential depth (i.e. Boltzmann constant 
k,B = 1). Geometric considerations for a three touching 
spheres configuration show that the choice of well-width 

5 = 0.5(v5 — 2\/3— 1) ~ 0.119 guarantees that each site 
is engaged at most in one bond. Hence, each particle can 
be form only up to two bonds and, correspondingly, the 
lowest energy per particle is — uq. 

The choice of a simple square-well interaction model 
to describe the association process between different par- 
ticles is particularly convenient from a theoretical point 
of view. It allows for a clear definition of bonding and a 
clear separation of the bond free energy in an energetic 
and entropic contributions, being unambiguously related 
to the depth of the well and to the bonding volume, re- 
spectively. 



in this class of theories, so that the modeling is not re- 
ally fully predictive (see discussion in Sect. Ill where this 
quantity is explicitly determined from Wertheim theory) . 
The unique aspect of W-TPT is that all the interaction 
parameters of this theory can be directly calculated from 
knowledge of the intermolecular potential and standard 
liquid state correlation functions, so that the theory is 
fully predictive. It has also been shown that W-TPT is 
formally equivalent to association-equilibrium models of 
self-assembly [20| . so that the Wertheim theory also of- 
fers the prospect of being able to improve the predictive 
character of these other theories if the theory itself can be 
validated as being reliable. The Wertheim theory itself is 
based on a formal perturbation theory [34 , l35l . l36j and 
there are naturally questions about the accuracy that can 
be expected from this theory. The present paper consid- 
ers a stringent test of Wertheim theory as a model of the 
thermodynamics of self-assembly by comparing precise 
numerical MC data for the thermodynamic properties of 
our model associating fluid to the analytic predictions 
of the Wertheim theory where there are no free parame- 
ters in the comparison. Notably, many of the properties 
that we consider have never been considered before in 
Wertheim theory. 



II. TWO PATCHY SITES PARTICLE MODEL 

We focus on a system of hard-sphere (HS) particles (of 
diameter tr, the unit of length) whose surface is decorated 
by M = 2 identical sites oppositely located (see Fig. [T]). 

The interaction V(l,2) between particles 1 and 2 is 



7(1,2) = V HS (r 12 ) + E E M r S) 



(1) 



i=l,2 .7=1,2 



where Vrs is the hard-sphere potential, Vw{x) is a 
square-well interaction (of depth — uq for x < 6, other- 
wise) and ri2 and r 1 ^ are respectively the vectors join- 
ing the particle-particle centers and the site-site (on dif- 
ferent particles) locations. Temperature is measured in 



III. WERTHEIM THEORY 

The first-order Wertheim thermodynamic perturba- 
tion theory [34], HH, H3 provides an expression for the 
free energy of associating liquids. The Helmholtz free 
energy is written as a sum of the HS reference free en- 
ergy Ahs plus a bond contribution A{, nd, which derives 
by a summation over certain classes of relevant graphs 
in the Mayer expansion (37} ■ In the sum, closed loops 
graphs are neglected. The fundamental assumption of 
W-TPT is that the conditions of steric incompatibilities 
are satisfied: (i) no site can be engaged in more than 
one bond; (ii) no pair of particles can be double bonded. 
These steric incompatibilities are satisfied in the present 
model thanks to the location of the two sites and the cho- 
sen 5 value. In the formulation of Ref. [38(, for particles 
with two identical bonding sites, 



N 



= 2\nX - X + 1 



(2) 



Here f3 = 1/UbT and X is the fraction of sites that are 
not bonded. X is calculated from the mass-action equa- 
tion 



X 



1 



1 + 2pXA 



(3) 



where p = N/V is the particle number density and A is 
defined by 



A = 4tt / 9H S (ri2)(f(l2)) u 



r\ 2 dri2 



(4) 



Here g#s(12) is the reference HS fluid pair corre- 
lation function, the Mayer /-function is /(12) = 
exp(-V w {Y\\)/kBT)-l, and (/(12)) Wlia)2 [39[ represents 
an angle-average over all orientations of particles 1 and 2 
at fixed relative distance ri 2 . Since all bonding sites are 
identical (same depth and width of the square-well in- 
teraction), A refers to a single site-site interaction. The 
number of attractive sites on each particle is encoded in 
the factor 2 in front of A in Eq. H In the W-TPT, the 
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resulting free energy is insensitive to the location of the 
attractive sites, i.e. to the bonding geometry of the par- 
ticle. Note that the angle averaged Mayer /-function co- 
incides with the bonding interaction contribution to the 
virial coefficient. At low T (i.e. {3u» 1) the hard-core 
contribution to the virial becomes negligible as compared 
to the bonding component. In this limit, it is also possible 
to assume exp(/3uo) -1» exp(/?uo), so that the averaged 
Mayer /-function can be approximated with the virial as 
well as with the integral of the Boltzmann factor over the 
bond volume [13]. 

For a site-site square- well interaction, the Mayer func- 
tion can be calculated as [39| 



(/(12)) UllWa = [exp(/?uo)-l]S(r) 



where 



Sir) 



(5 + <7-r) 2 (25-a + r) 
6a 2 r 



(5) 



(6) 



is the fraction of solid angle available to bonding when 
two particles are located at relative center-to-center dis- 
tance r. Thus the evaluation of A requires only an ex- 
pression for gHs{ T i2) in the range where bonding occurs 
(a < r < a + 5). We have used the linear approxima- 
tion mi 



9Hs{r) 



1-0.50 9 0(1 
(1-0)3 "2(T= 



(7) 



(where = ^a 3 ' p) which provides the correct Carnahan- 
Starling [I^] value at contact. This gives 



V b (e 



0u o 



1) 



5 (3 + 86 + 35 2 ) 
2 (15 + 45) 



(1 - 0)3 

3 (125 + 5S 2 ) 
2 (15 + 45) 



(8) 



where we have defined the spherically averaged bonding 
volume V b = 4tt f° +S S(r)r 2 dr = tt<5 4 (15 + 4<5)/30. For 
the specific value of S studied here, Vb = 0.000332285 a 3 . 
At low 0, gHs(r) tends to the ideal gas limit value 
9Hs( r ) ~ 1. In this limit 



A = V h (e (iu ° - 1) 



(9) 



Eq. Q can be easily solved, providing the T and p 
dependence of X as 



X = 



1 + Vl + 8pA' 



(10) 



which has a low T limit X as \J2pA. 

In the more transparent chemical equilibrium fori 
Eq. (fTU|) can be written as 



1 - X 



2pA = pK b . 



(11) 



The last expression shows that, within Wertheim the- 
ory, bonding can be seen as a chemical reaction between 
two unreacted sites forming a bonded pair. In this lan- 
guage the quantity Kb = 2A plays the role of equilib- 
rium constant (in unit of inverse concentration) . Writing 
pKb = exp{— f3(AUb — TASb)} (introducing the energy 
and entropy change in the bond process), it is possible to 
provide precise expressions for AUb and ASb within the 
Wertheim theory. Specifically, when e^ u ° >> 1 (a very 
minor approximation since aggregation requires T « uq 
to be effective) it is possible to identify 



AS h = In 



8np 



a+S 



AU b = -u (12) 
g HS (r)S{r)r 2 dr (13) 



(note that ASb is in dimensionless entropy units, since 
Boltzmann's constant is taken to equal 1.). If gns ~ 
1, ASb — ln(2NVb/V). Hence the change in energy is 
given by the bond energy, while the change in entropy is 
essentially provided by the logarithm of the ratio between 
the bonding volume and the volume per site (V/2N). 



IV. CLUSTER SIZE DISTRIBUTIONS AND 
ASSOCIATION PROPERTIES 

It is interesting to discuss the prediction of the 
Wertheim th eory in term of clusters of physically bonded 
particles [43|, [4J] . In the case of square- well interactions 
(differently from a continuous potential) we can define a 
bond between two particles unambiguously. Evidently, 
when there is a bond the interaction energy equals - uo- 

To make the discussion more transparent, we can de- 
fine pb = 1 — X as the probability that an arbitrary site 
is bonded. It is thus easy to convince oneself that the 
number density of monomers is p\ = p(l — pb) 2 — pX 2 , 
since both sites must be unbonded 0. Similarly a chain 
of I particles has a number density p\ equal to 



Pl = P {i-Pb?pt 1) = P x\i-xf-v 



(14) 



since one site of the first and one of the last particle 
in the chain must be unbonded and I — 1 bonds link the 
I particles. 

Once the cluster size distribution of chains is known, 
it is possible to calculate the average chain length L as 
the ratio between the total number density and the num- 
ber density of chains in the system, (i.e., as the ratio 
between the first (I 1 ) and the zero (1°) moments of the 
pi distribution) 



i 

x 



(15) 



where we have substituted, by summing the geometric 
series over all chain lengths, X^z=i I Pi = P an( l E;=i Pi — 



4 



pX. Thus, using Eq. (JTOJI, 

_ 1 + VI + 8pA 



(16) 



At low T, L w ^2pA and hence L grows in density as 
y/p (if the density dependence of A can be neglected [451 ]) 
and in T as L ~ exp(/3uo/2). 

The potential energy of the system coincides with the 
number of bonds (times — Uo). Hence, the energy per 
particle E/N is 



E/N = -u 



-uq(1 — X) = -uopb (17) 



J2Zi 1 pi 

The same result is of course obtained calculating E/N as 

d(f3A b ^ d /N) where ^ AbQnd / N ) is given by Eq . Thc 

energy approaches its ground state value (E gs /N = —uq) 
as 



AT 



m X w u (2pA) 2 



(18) 



i.e. with an Arrhenius law with activation energy Uo/2 
in the low T limit, a signature of independent bonding 
sites 0,113, S|. 

From Eq. (TIT)) it is possible to calculate the (constant 
volume) specific heat Cy as 



C V = C 



X 2 e f3uo p 
VI + 8pAT 2 



(19) 



where C — 87rwn 



9Hs{ r )S{r)r 2 dr. In thc low T 




FIG. 2: W-TPT predictions for the T-dependence of the spe- 
cific heat Cv for different values of densities p (Eq. I19[) . The 
inset shows the value of the specific heat at the maximum 



nite T (see Fig. [5]), which define a lines of specific heat 



extrema in the T — p plane. The location of the max- 
imum in the specific heat has been utilized to estimate 
the polymerization temperature [13, 0, d, H3, 0, 

Within the theory, it is also possible to evaluate the 
extent of polymerization <I>, defined as the fraction of 
particles connected in chains (chain length larger than 
one), i.e. 



$ = 



IX2 lPl 

J2Zi 1 pi 



= 1 - 



pi 
p 



(20) 



where we have used Eq. [TH <I> plays the role of order 
parameter for the polymerization transition. The density 
and temperature dependence of $ are shown in Fig. [3] 
The cross-over from the monomeric state at high T ($ « 

0) to thc polymeric thermodynamic state at low T ($ « 

1) takes place in a progressively smaller T- window on 
decreasing p 
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FIG. 3: Extent of polymerization $ vs. temperature T (a) 
and density p (b). Symbols are GC simulation data. 

Indeed, for each p, a transition temperature T$ can 
be defined as the inflection point of $ as function of 
T [13, [H, [H, H H3- The locus T$(p) provides an es- 
timate of the polymerization line in the phase diagram. 
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The location of the inflection point of the energy E and 
of $ are different, since dE/dT ~ dX/dT (Eq.Q7|), while 
d$/dT XdX/dT (Eq.HJJ). In other models of equi- 
librium polymerization, incorporating thermal activation 
or chemical initiation 7$ and Tc™--* coincide (49l [5l| . 

Another estimate of the transition line of this rounded 
thermodynamic transition can be defined as the locus in 
the T — p plane at which $ = 0.5, i.e. half of the parti- 
cles are in chain form (the analog of the critical micelle 
concentration [12]), i.e. /3i(T 1 / 2 ,p) = p/2. The corre- 
sponding temperature T x / 2 is then given by the solution 
of the equation X = \/0.5, or equivalently 

pA = 1 - (21) 

In the present model, the <& = 0.5 locus is a line at con- 
stant pi,, corresponding to a constant value of the product 

p [exp(/3u ) - !]• 

To provide a global view of the polymerization transi- 
tion, we show in Fig. 0] the Wertheim theory predictions 
for the specific heat maximum and polymerization transi- 
tion lines. Curves becomes progressively more and more 
similar on cooling. As clearly shown by the simple ex- 
pression for T1/2 (Eq. I21[) the quantity pA is constant, 
which implies that lnp ~ l/T. This behavior is also ap- 
proximatively found for the loci defined by the inflection 
point of <E> and E, as shown in Fig. [TJ 
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FIG. 4: Specific heat maxima and polymerization transition 
lines, as predicted by the Wertheim theory, in the T — p plane. 
The inset shows the linear scale. Observe the similarity of the 
data in Fig. 2-4 to that of Fig. 4 - 6 in Ref. [2^] correspond- 
ing to equilibrium polymerization in the Stockmayer fluid, a 
Lennard- Jones particle with a superimposed point dipole. 

It is interesting to evaluate the value of average degree 
of polymerization L along the polymerization transition 
line. This information helps estimating the polymeriza- 
tion transition itself, but it is also relevant for the re- 
cently proposed analogies between polymerizing systems 
and glass-forming liquids [Hj]. As shown in Fig. [5] the 
transition takes place for L « 2, consistent with previous 
findings (compare with the inset of Fig. 7 in Ref. for 
the Stockmayer fluid). The fact that, at T^, L ss 2 can 



provide a way to locate the transition temperature when 
experimental (or numerical) data are noisy [32| . 
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FIG. 5: Value of the average chain length along the polymer- 
ization lines T*(p) and Tc™ ax (p), according to the Wertheim 
theory. 

A final relevant consideration concerns the expression 
for the bonding free energy. Under the assumption of an 
ideal state of chains, the system free energy A can be 
written as 

oo 

/3A/V = Y / pi\hip l -l + {l-l)hiK b ] (22) 
i=i 

which accounts for the translational entropy of the clus- 
ters as well as for the bond free energy In Kb which 
is assumed to be linear in the number of bonds. At 
high T only monomers are present, pi — pSn and hence 
(3A/V = p[ln p — 1], The bonding part of the free energy 
i3Ab on d can thus be evaluated as difference of (3 A and the 
corresponding high T limit, i.e. 

oo 

pAbond/V = Y^Pi [In Pi - 1 + (i - 1) m^h] - ppnp - 1] 

i=i 

(23) 

Inserting the Wertheim cluster size distribution (Eq. [Tl)l . 
summing over all cluster size and using the definition 
Kb = 2 A (Eq. ITT]) . one exactly recovers the Wertheim 
bonding free energy (Eq. [5]). The Wertheim theory 
can thus be considered as a mean-field theory of chain 
associations, with the hard-sphere free energy as its 
high-T limit. Differently from other mean-field ap- 
proaches [13, EE El) H3 j the theory provides also a well- 
defined prescription for calculating the equilibrium con- 
stant, which partially account for the structure of the 
reference hard-sphere fluid via the gnsif) contribution 
in Kb- It is this special feature that makes the theory 
fully predictive. 

V. MONTE CARLO SIMULATION 

We have performed numerical simulations of the model 
in the grand-canonical (GC) ensemble [55j for several val- 
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ues of T and of activities to evaluate the structural prop- 
erties of the system as a function of density and temper- 
ature. We have performed two types of grand canonical 
simulations: particle and chain insertion/removal. 

The first method (particle) is a classical GC MC simu- 
lation where monomers are individually added to or elim- 
inated from the system with insertion and removal prob- 
abilities given by 

zV 

Pinsertion = mm(l, -e^ 25 ) (24) 

N e^ AB 

Premoval = min(l, — ), (25) 

V z 

where N is the number of monomers, AE is the change 
in the system energy upon insertion (or removal) and z 
is the chosen activity. We have simulated for about two 
million MC steps, where a MC step has been defined as 
50000 attempts to move a particle and 100 attempts to 
insert or delete a particle. A move is defined as a displace- 
ment in each direction of a random quantity distributed 
uniformly between ± 0.05 a and a rotation around a ran- 
dom axis of random angle distributed uniformly between 
±0.1 radiant. The box size has been fixed to 50 a. With 
this type of simulations we have studied three tempera- 
tures (T = 0.08, 0.09 and 0.1) and several densities rang- 
ing from 0.001 up to 0.2 (corresponding to number N of 
particles ranging from N w 1000 to N « 20000. 

The second method is a grand-canonical simulations 
where we had and remove chains of particles (see for ex- 
ample Ref. (HI). The insertion and removal probabilities 
for a chain of length I are: 

z l e -(l-l)f3f b y 

JWtt* = mi»(l, ^ e-' ) (26) 

Premoval = min(l, — ^^(j^)^ ) (27) 

where Ni is the number of chains of size I, AE is the 
change in energy and (3f is the In of the integral of the 
Boltzmann factor over the bond volume, i.e., 

Pfb = 0u o - H2V b ] (28) 

The activity of a chain of / particles is thus written as 
z% = z l e~U~ 1 )Pf b , or equivalently in term of chemical po- 
tential Hi = l/j, — (I — l)/b, where p, is the monomer chem- 
ical potential. Using the chain insertion algorithm, we 
have followed the system for about 10 5 MC steps, where 
a MC step has now been defined as 50000 attempts to 
move a particle (as in the previous type) and 100 at- 
tempts to insert or delete a chain of randomly selected 
length. The geometry of the chain to be inserted is also 
randomly selected. The box size was varied from 50 to 
400 a, according to density and temperature, to guar- 
antee that the longest chain in the system was always 
shorter than the box size. At the lowest T, TV « 50000. 
Using chain moves we have been able to equilibrate and 




FIG. 6: Snapshots of a fraction of the simulated system at 
p « 0.0035 and T = 0.10 and T = 0.055. At this density, the 
polymerization transition is located at T® ps 0.08. The length 
of the shown box edge is about 30a. 



study densities ranging from 10~ 6 up to 0.2 for four low 
temperatures (T = 0.05, 0.055, 0.6 and 0.7) where chains 
of length up to 400 monomers are observed. We have 
also studied the same three temperatures (T = 0.08, 0.09 
and 0.1) examined with the particle insertion/removal 
method to compare the two MC approaches. 

As a result of the long simulations performed (about 
three months of CPU time for each state point) resulting 
average quantities calculated from the MC data (L, $, 
E/N) are affected by less than 3 % relative error. Chain 
length distributions pi (whose signal covers up to six or- 
der of magnitude) are affected by an error proportional 
to \log(pi) which progressively increases on decreasing pi, 
reaching 70% at the smallest reported pi values. 



VI. SIMULATION RESULTS 

First, we begin by showing in Fig. [^representative par- 
ticle configurations both above and below the polymer- 
ization transition temperature, T^. Evidently, the parti- 
cles are dispersed as a gas of monomers and as a gas of 
chains above and below this characteristic temperature. 
The low-T configurations is composed by semi-flexible 
chains, with no rings. Indeed, the short interaction 
range 5 introduces a significant stiffness in the chain and 
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chain length 1 



FIG. 7: Mean squared end-to-end distance < > for iso- 
lated chains of length I built according to the model studied 
in this article. 



a persistence length extending over several monomers. 
To quantify the linearity of the chains for the present 
model we show in Fig. [7] the chain end-to-end squared 
distance < J?f > for isolated chains of chain length up to 
I ~ O(10 3 ). Single chains are generated by progressively 
adding monomers to a pre-existing chain in a bonding 
configurations, after checking the possible overlap with 
all pre-existing monomers. Since the bond interaction is a 
well, all points in the bond-volume have the same a-priori 
probability. As shown in Fig. [7j the end-to-end chain dis- 
tance scales as a power-law (< B? e >~ l 2v ) both at small 
and large I values, with a crossing between the two behav- 
iors around I ~ O(10). At small I, the chain is persistent 
in form and thus is rod-like. For larger I, the best-fit with 
a power law suggests an apparent exponent 2v m 1.1, 
that is expected to evolve — for very lon g ch ains — to- 
ward the self-avoiding value 2v w 1.18 (56l.l57||. We recall 
that 2v = 6/5 in the Flory mean field prediction [58| and 
2v — 1 in the simple random walk model. 

To provide evidence that the chain GC MC simulation 
provides the correct sampling of the configurations (and 
hence that the activity of the chain of length I is correctly 
assigned) we compare in Fig. [8] the chain length densities 
pi calculated with the two methods at T = 0.08. The dis- 
tributions calculated with the two different methods are 
identical. Similar agreement is also found at T = 0.09 
and T = 0.1. This strengthens the possibility of using 
the chain MC method, which does not significantly suf- 
fer from the slow equilibration process associated to the 
increase of the Boltzmann factor exp (/3uo) on cooling. 
All the following data are based on chain GC-MC simu- 
lations. 

We next compare the chain length distributions calcu- 
lated using the chain MC method with the predictions of 
the Wertheim theory. We compare the simulation data 
with two different approximation: in the first one we 
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FIG. 8: Comparison of the chain length distributions pi 
at T — 0.08 obtained independently with the monomer 
(full symbols) and chain (lines connecting open symbols) 
grand-canonical simulations for several values of the monomer 
activity z. From left to right, the activity values are: 
10 _3 ,1.5 x 10~ 3 ,2.4 x 10~ 3 ,3 x 10 _3 ,5 x 10 _3 ,6 x 10 -3 . 
Similar agreement is also found at the other temperatures. 



choose the ideal gas as reference state, i.e. we approxi- 
mate the reference radial distribution function with one. 
In the more realistic approximation, we use the small r 
expansion of the hard-sphere radial distribution function 
(see Eq. [7]). Comparison between simulation data and 
theoretical predictions (note there are no fitting parame- 
ters) is reported in Fig. [5] for two different temperatures. 
At low densities (sampled at low T) the approximation 
gjjs ~ 1 is already sufficient to properly describe p;. At 
higher densities (sampled at higher T), the full theory is 
requested to satisfactory predict the chain length distri- 
butions. 

Fig. HOf a) compares the Wertheim theory predictions 
for the average chain length with the corresponding sim- 
ulation results. In the entire investigated p and T range, 
the Wertheim theory provides an accurate description 
of the equilibrium polymerization process. The limiting 
growth law in y/p is clearly visible at the lowest tem- 
peratures. At the highest temperatures, it is possible to 
access the region of larger densities (p ~ 0.1) where the 
presence of other chains can not be neglected any longer 
and A becomes p dependent. In this limit, the growth 
low L ~ yfp is no longer obeyed [45] . As a further check 
on Wertheim theory, we collapse all the L data to the 
universal functional form predicted by Wertheim theory 
using the scaling variable 2Ap (Fig. flOl b)) 

As an ulterior confirmation of the predictive capabili- 
ties of Wertheim theory, we report a comparison between 
simulations and theory for the density dependent of the 
extent of polymerization <t> (Fig. [3J and for the energy 
per particle (Fig. ITT|) . Both figures clearly shows a excel- 
lent agreement between the simulated and the predicted 
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FIG. 9: Chain length densities pi for several activity values at 
two different temperatures, T = 0.06 and T = 0.09. Symbols 
are simulation data. Lines are Wertheim theory predictions: 
dashed lines assume gus ~ 1 (Eq. [9] for A), while full lines 
are based on the full radial dependence of qhs (Eq. [5] for 
A). At the lowest T (and average densities), the ideal gas 
approximation is sufficient to model the Monte Carlo data. 



p dependence of <& and E at all T investigated. 

As a test of the validity of the approach of the poly- 
merizing system as an ideal gas of equilibrium chains, we 
compare the monomer activity z and the monomer den- 
sity pi in Fig. [T2J Indeedthe activity of a cluster of size I 
coincides with pi , which is consistent with ideal gas scal- 
ing. In particular, the activity of the single particle (the 
input in the MC grand-canonical simulation) can be com- 
pared with the resulting density of monomers p\ . Data in 
Fig. [T2"l shows that the ideal-gas law is well obeyed at low 
T and p, confirming that, in the investigated range, the 
system can be visualized as an ideal gas mixture of chains 
of different lengths, distributed according to Eq. [T4"l 

The existence of a large T-p window where an ideal 
mixture of chains provides a satisfactory representation 
of the system suggests that, in this window, correlations 
between different chains can be neglected. In this limit, 
the structure of the system should be provided by the 
structure of a single chain, weighted by the appropriate 
chain length distribution. Specifically, we have 
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FIG. 10: Average chain length as a function of the density for 
all studied temperatures, (a): Lines are the Wertheim theory 
predictions: dashed lines assume gas = 1 (Eq.|9]for A), while 
full lines are based on the full radial dependence of <?£rs(Eq.|8] 
for A), (b): Scaled representation of L vs. 2pA. Symbols 
are simulation data. The line is the function L(x) = 1 +^^+ 2x 
(See Eq. [IB) 



S(q) = 



1 E> 



iq-(r 



(29) 



where r» is the coordinate of particle i, the sum runs over 
all N particles in the system, and the average (...) is over 
equilibrium configurations. In the ideal gas limit, cor- 
relations between different chains can be neglected and 
S(q) can be formally written as, 



S(q) 



EZi pi 1 



(30) 



where Si(q) is the structure factor (form factor) of a chain 
of length I: 



Si(q) 



i 

^ e -if-(fi 



(31) 
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FIG. 11: Energy per particle in unit of Mo as a function of the 
density for all studied temperatures. Symbols are simulations 
data while lines are the Wertheim theory predictions: dashed 
lines assume gns = 1 (Eq. [§]for A), while full lines are based 
on the full radial dependence of gns (Eq. [8] for A) 
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where < l m >= J^i l m Pl denotes the m moment of the 
cluster size distribution pi . Fig. [T3] shows a comparison 
between the S(q) calculated in the simulation and the 
theoretical S(q) evaluated according to Eq. I3TH and I3"2l at a 
low T, where the ideal gas approximation is valid. Devia- 
tions are only observed at the highest density, suggesting 
that the ideal gas of chains is a good representation of 
the structure of the system, in agreement with the equiv- 
alence between activity and monomer density shown in 
Fig. [TJl This observation is particularly relevant, since 
it suggests that (in the appropriate T-p window) also a 
description of the dynamics of the model based on the 
assumption of independent chains can be attempted. 
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FIG. 13: Comparison of the structure factor calculated from 
the simulated configurations and from Eq. I30ll32l at T=0.055 
for different values of the activity. 



FIG. 12: Relation between the activity z and the monomer 
density p\ for all investigated temperatures. Line indicates 
the ideal gas law z = pi . 



Since the persistence length of the chains is rs 10 — 20 
particles (see Fig. [7J) , one can assume, as a first approx- 
imation, that in the investigated T-p region chains are 
linear. When this is the case, averaging over all possible 
orientation of the chain gives: 



Sl(q) 



l + i^2(/-i) 



3=1 



sin(jqa) 



The small q expansion of Si(q) is 

l(l 2 1) 



Si(q) « I 



:-!(» 



Correspondingly, S(q) behaves at small q as 



S(q) 



<l 2 > 
< I > 



1 - 



q 2 a 2 
36 



< / 4 > - < I 2 > 
<l 2 > 



(32) 



(33) 



(34) 



VII. CONCLUSIONS 

In our pursuit of a fully predictive molecularly-based 
theory of self-assembly in terms of molecular parameters 
and liquid state correlation functions, we have considered 
a direct comparison of a liquid in which fluid particles 
have sticky spots on their polar regions to the predic- 
tions of the Wertheim theory for the relevant properties 
governing the self-assembly thermodynamics. In the in- 
vestigated region of temperatures and densities (which 
spans the polymerization transition region), the predic- 
tions of Wertheim theory describe simulation data re- 
markably well (without the use of any free parameters 
in these comparisons). This success means that we can 
be confident in pursuing more complicated types of self- 
assembly based on the foundation of Wertheim theory. 
For example, it is possible to extend the Wertheim the- 
ory analysis in a direct way to describe the self-assembly 
of branched chains having multifunctional rather than 
dipolar symmetry interactions as in the present paper 
and to compare these results in a parameter free fash- 
ion to corresponding MC simulations for the these mul- 
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tifunctional interaction particles [591 ]. Recent work has 
shown that the Wertheim theory describes the critical 
properties of these mutifunctional interaction particles 
rather well [60l ]. According to this recent study, liquid 
phases of vanishing density can be generated once small 
fraction of polyfunctional particles are added to chain- 
forming models like the one studied here. With the new 
generation of non-spherical sticky colloids, it should be 
possible to realize "empty liquids" [(5(| and observe equi- 
librium gelation [46l [6l|. i.e. approach dynamical arrest 
under equilibrium conditions. 

The Werhtheim theory has also been applied success- 
fully to description of molecular associated liquids [62l.l63l. 
[rjij ] and to the thermodynamics of hard sphere polymer 
chains with short range attractive interactions [65l [66[ . 
Thus, the theory could be adapted to describing mutu- 
ally associating polymers and the formation of thermally 
reversible gels in these fluids upon cooling. 

In summary, the Wertheim theory provides a promis- 
ing framework for treating the thermodynamics of a wide 
range of self-assembling systems. The development of 
this theory and its validation by simulation and measure- 
ment should provide valuable tools in the practical devel- 
opment of self-assembly as a practical means of synthetic 
manufacturing. This theory also offers the prospect of 
improving the existing equilibrium association theories 
that are largely based on a lattice fluid model framework. 
This could allow progress to be made more rapidly, since 
this type of computation often offers computational ad- 
vantages and because many problems such as chemically 
initiated chain branching [67| and thermally activated as- 
sembly processes have already been considered by lattice 
approaches |5Q. |68|. 

The problem of estimating the entropy of association 
in real self-assembling molecular and particle systems in 



solution is a difficult problem that has been addressed by 
many authors previously (Ref. [6^| and refs. therein). It 
would clearly be interesting to extend the present work 
to determine how well the Wertheim theory could pre- 
dict entropies of association for self-assembly processes 
that occur in a solvent rather than in the gas phase. The 
most interesting solvent in this connection, water, is a 
particular challenge since water itself can be considered 
an associating fluid, so that we are confronted with the 
problem of how the water association couples to the par- 
ticle self-assembly. The problem of understanding the 
common tendency of particle self-assembly in aqueous 
solutions to occur upon heating requires particular in- 
vestigation. In the future, we look forward to exploring 
these more complex mixtures of associating fluids, which 
are so prevalent in real biological systems and in a ma- 
terials processing context. 

As a final comment, we note that numerical work on 
this class of simple models (playing with the particle in- 
teraction symmetries) can help understanding more com- 
plicated ordered structures (as for example sheet-like, 
nanotube and closed nanoshell structures), as recently 
found when particles have multipole interaction poten- 
tials [1, [t| 0j], EH- We also note that the results dis- 
cussed here apply to the growing field of functionalized 
colloidal particles, colloidal particles with specifically de- 
signed shapes and interaction sites 0, 0, d, 0, [H, [7(| • 
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